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#14  ABSTRACT 

We  developed  a  novel  approach  to  extend  the  particle  level  set  method  to  the  simulation 
of  as  many  regions  as  desired.  The  various  regions  can  be  liquids  or  gases  of  any  type 
with  differing  viscosities,  densities,  viscoelastic  properties,  etc.  We  also  proposed 
techniques  for  simulating  interactions  between  materials,  whether  it  be  simple  surface 
tension  forces  or  more  complex  chemical  reactions  with  one  material  converting  to 
another  or  two  materials  combining  to  form  a  third.  When  discretizing  the  underlying 
Navier-Stokes  equations  for  multiphase  flow,  an  additional  difficulty  occurs  since 
discretization  stencils  cross  region  boundaries  naively  combining  non-smooth  or  even 
discontinuous  data.  Recently,  we  developed  a  new  coding  paradigm  that  allows  one  to 
incorporate  physical  jump  conditions  in  data  "on  the  fly,"  which  is  significantly  more 
efficient  for  multiple  regions  especially  at  triple  points  or  near  boundaries  with  solids. 
This  removes  the  need  for  any  algorithm  changes  that  might  reduce  the  accuracy  of  the 
scheme,  and  moreover  even  removes  the  need  for  changes  to  the  code  itself.  Besides  this 
work  we  have  also  addressed  scalability  including  methods  on  octree  and  Run  Length 
Encoded  (RLE)  data  structure,  as  well  parallel  implementation  such  as  MPI.  Other  work 
includes  work  on  fracture. 
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The  particle  level  set  method  [EFFM,EMF]  was  developed  under  prior  ONR  funding,  and 
has  proven  successful  for  the  numerical  simulation  of  two  separate  regions.  In  the  case  of 
free  surfaces,  the  Dirichlet  boundary  conditions  imposed  on  the  water/air  boundary  were 
increased  to  second  order  accuracy  in  [ENGF]  using  the  technique  proposed  in  [GFCK], 
Later,  [ELF]  showed  that  the  level  set  advection  could  be  handled  with  a  simple  first 
order  accurate  semi-Lagrangian  method  as  long  as  the  particles  were  adverted  with 
second  order  accuracy.  This  allowed  for  the  extension  of  the  method  to  more  complex 
phenomena,  such  as  coupling  to  thin  deformable  objects  as  in  [GSLF].  This  plays  a  key 
role  in  a  recent  MURI  in  which  Farhat  and  Fedkiw  compose  the  Stanford  team.  This 
work  on  coupling  two  phase  flows  to  both  rigid  and  deformable  solids  leveraged  our 
recent  algorithms  for  rigid  body  simulation  [GBF],  and  finite  element  simulation  of 
deformable  volumetric  bodies  [ITF]  and  thin  shells  [BFA].  We  even  developed  numerical 
techniques  that  allowed  a  deformable  solid  modeled  with  a  Lagrangian  finite  element 
mesh  to  phase  change  into  a  liquid  modeled  with  an  Eulerian  background  mesh. 
However,  all  of  these  algorithms  were  still  limited  to  at  most  two  liquid  phases. 

Thus,  most  recently  we  developed  a  novel  approach  to  extend  the  particle  level  set 
method  to  the  simulation  of  as  many  regions  as  desired  [LSSF].  The  various  regions  can 
be  liquids  or  gases  of  any  type  with  differing  viscosities,  densities,  viscoelastic 
properties,  etc.  We  also  proposed  techniques  for  simulating  interactions  between 
materials,  whether  it  be  simple  surface  tension  forces  or  more  complex  chemical 
reactions  with  one  material  converting  to  another  or  two  materials  combining  to  form  a 
third.  We  use  a  separate  particle  level  set  method  for  each  region,  and  propose  a  novel 
projection  algorithm  that  decodes  the  resulting  vector  of  level  set  values  providing  a 
"dictionary"  that  translates  between  them  and  the  standard  single-valued  level  set 
representation. 

When  discretizing  the  underlying  Navier-Stokes  equations  for  multiphase  flow,  an 
additional  difficulty  occurs  since  discretization  stencils  cross  region  boundaries  naively 
combining  non-smooth  or  even  discontinuous  data.  Under  previous  ONR  funding,  we 
introduced  the  ghost  fluid  method  [FAMO]  some  years  ago  in  order  to  treat  this  problem 
for  compressible  flow.  Sine  then,  it  has  been  extended  to  incompressible  flows 
[LFK,KFL]  and  combustion  [NFK,NFJ]  and  addressed  by  many  other  authors.  Although 
the  original  idea  in  [FAMO]  was  to  create  ghost  cells,  many  authors  have  explored 
various  algorithm  or  implementation  ideas  to  avoid  the  ghost  cells  with  the  aim  of 
achieving  identical  or  at  least  similar  results.  In  [LSSF],  we  developed  a  new  coding 
paradigm  that  allows  one  to  incorporate  physical  jump  conditions  in  data  "on  the  fly," 
which  is  significantly  more  efficient  for  multiple  regions  especially  at  triple  points  or 
near  boundaries  with  solids.  This  removes  the  need  for  any  algorithm  changes  that  might 
reduce  the  accuracy  of  the  scheme,  and  moreover  even  removes  the  need  for  changes  to 


1 


the  code  itself.  A  smart  discretization  class  is  used  to  encode  both  the  region  one  is 
working  in  as  well  as  the  required  jump  conditions,  and  the  programmer  simply  writes 
the  code  as  if  it  were  for  a  single  region  relying  on  the  smart  discretization  class  to  create 
the  correct  code. 

Computational  resources  are  almost  always  limited,  and  this  problem  cannot  typically  be 
addressed  by  getting  more  machines  as  the  algorithms  do  not  always  scale  to  multiple 
processors.  Thus,  there  is  a  variety  or  research  into  way  of  improving  the  efficiency  of 
the  algorithms  themselves  such  as  adaptive  mesh  technology.  In  fact,  under  ONR 
funding,  we  proposed  first  algorithms  for  multiphase  flow  simulation  on  an  octree  data 
structure  in  [LGF].  Later  in  [LGO],  we  both  extended  our  discretization  to  second  order 
accuracy  as  well  as  proposing  a  more  efficient  technique  that  uses  a  top  level  Cartesian 
grid  placing  an  individual  octree  in  each  Cartesian  cell.  Most  recently  we  have  examined 
methods  for  discretizing  large  bodies  of  water  where  researchers  typically  use  either  the 
deep  water  or  shallow  water  approximations.  We  instead  wanted  to  use  a  fully  three 
dimensional  technique  to  capture  wave  overturning  and  related  effects,  but  obviously  had 
to  use  an  adaptive  mesh  because  of  computational  limitations.  Instead  of  using  the  octree 
or  a  patch  based  AMR  approach,  we  designed  a  new  method  based  on  Run  Length 
Encoding  (RLE)  from  image  processing.  RLE  takes  an  image  of  pixels,  and  compresses 
them  by  storing  row  subsets  with  similar  properties  in  a  smaller  data  structure.  Applying 
the  same  idea  to  the  incompressible  Navier-Stokes  equations,  we  start  with  a  uniform 
Cartesian  grid  and  merge  any  two  or  more  cells  in  order  to  make  tall  cells.  Note  that  we 
only  merge  cells  in  the  vertical  direction,  since  we  will  be  treating  pressure  linearly  in 
each  of  these  cells,  and  a  linear  pressure  profile  is  more  valid  in  water  than  a  linear 
horizontal  one.  In  fact,  if  every  column  of  water  were  merged  into  one  cell,  we  would 
obtain  a  two-dimensional  height  field  approximation  similar  to  the  shallow  water 
equations.  Thus  we  typically  merge  parts  of  columns  under  water  while  leave  fully  three- 
dimensional  approximations  near  the  surface  of  the  water  and  near  any  objects  (which 
may  be  submersed).  Starting  with  the  standard  MAC  grid  discretization  of  the 
incompressible  Navier-Stokes  equations,  [IGLF]  develops  the  natural  discretizations 
obtained  by  merging  cells  in  this  vertical  fashion  assuming  linear  pressure  profiles.  A 
current  limitation  of  the  merged  tall  cells  is  that  they  cannot  adequately  resolve  vorticity. 
However,  we  have  also  recently  developed  a  hybrid  simulation  framework  that  leverages 
a  vortex  particle  approach  on  a  background  grid  [SRF],  This  method  can  also  be 
exteneded  to  properly  carry  vorticity  through  the  merged  tall  cells. 

In  addition  to  adaptive  mesh  strategies,  we  have  been  working  on  ways  to  make 
algorithms  more  scalable.  Recently,  our  group  was  awarded  a  DURIP  which  was  used  to 
purchase  an  80  processor  SUN  with  600G  of  RAM.  Since  obtaining  the  cluster,  we  have 
worked  to  parallelize  our  code  base  using  MPI  and  currently  have  the  full  fluid 
simulation  system  parallelized.  Also,  both  the  PI  and  two  graduate  students  involved  in 
this  work  received  screen  credits  on  the  film  Poseidon  -  which  has  been  nominated  for  an 
Oscar  for  special  effects  work.  For  that  film,  Industrial  Light  +  Magic  made  used  of  our 
new  MPI  code  base  to  simulate  all  of  the  ocean  scenes.  In  that  film,  the  entire  ship  and 
ocean  were  completely  computer  generated  using  our  PhysBAM  software.  Of  course,  it  is 
not  enough  to  merely  write  the  MPI  code,  but  rather  one  needs  to  find  the  bottlenecks  and 
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improve  the  implementation  -  or  even  create  better  algorithms.  Thus,  we  have  been 
collaborating  with  Intel  for  some  time  working  on  getting  algorithms  from  our  PhysBAM 
code  based  implemented  in  their  multicore  framework.  Two  students  from  our  group 
spend  one  day  at  Intel  on  a  weekly  basis.  Intel  is  pushing  for  a  suite  of  new  multicore 
benchmarks  and  our  code  base  will  be  part  of  this.  Currently,  we  our  working  on 
extending  our  finite  element  solids  simulations  to  the  MPI  frame  work  as  well. 

Previously,  we  have  done  quite  a  bit  of  work  on  the  finite  element  simulation  of  solids, 
some  of  which  is  mentioned  above.  The  early  work  focused  on  using  level  set  methods  to 
create  red-green  adaptive  tetrahedral  and  surface  triangle  meshes  for  subsequent  finite 
element  simulation  [MBTF,BTMF].  Our  main  interest  was  in  contact  and  collision,  as 
well  as  materials  that  undergo  large  deformation,  plasticity  and  fracture.  Thus,  a  main 
concern  was  that  the  finite  element  technique  would  be  robust  under  large  deformation 
and  we  devised  a  new  algorithm  in  [ITF]  that  not  only  allowed  tetrahedral  to  invert,  but 
that  provided  forces  that  would  work  to  un-invert  them.  When  materials  fracture,  it 
becomes  necessary  to  generate  a  new  mesh  for  the  finite  element  analysis.  Thus,  we 
proposed  a  new  virtual  node  algorithm  to  automatically  generate  this  mesh  in  [MBF]. 
This  method  has  a  similar  flavor  to  the  ghost  fluid  method  where  ghost  nodes  are  created 
to  model  discontinuities  in  the  material.  The  main  contribution  of  the  algorithm  is  that  it 
figures  out  exactly  how  to  most  efficiently  treat  all  the  degrees  of  freedom  by  creating 
ghost  nodes  as  parts  of  new  tetrahedral.  As  far  as  this  aspect  of  the  algorithm  is 
concerned,  it  turns  out  that  XFEM  is  a  subset  of  our  approach.  Once  these  new  degrees  of 
freedom  are  created,  one  can  proceed  in  almost  standard  fashion  with  a  robust  and 
efficient  finite  element  mesh.  Most  recently,  in  [BHTF],  we  have  extended  our  prior  work 
on  fracture  to  treat  truly  rigid  bodies  that  can  fracture  into  multiple  rigid  bodies.  Some  of 
the  key  aspects  of  the  approach  included  a  method  to  properly  handle  the  null  space  in  the 
rigid  body  finite  element  simulation  as  well  as  techniques  for  treating  the  subsequent 
contact  and  collision. 
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